
function paiOpt = GetTimingPrem(outEstim,pai0,numSim,pathLength)

model  = outEstim.model;
params = model.params;

% For anti-thetic sampling
rng(1,'twister');
shocksMat                     = zeros(size(model.eta,2),pathLength,numSim);
shocksMat(:,:,1:numSim/2)     = randn(size(model.eta,2),pathLength,numSim/2);
shocksMat(:,:,numSim/2+1:end) = -shocksMat(:,:,1:numSim/2);

% settings
setupTP.model  = model;
setupTP.params = model.params;
posVf          = find(strcmp(model.labely,'$vf_t$'));
setupTP.vf     = exp(model.g0(posVf));

% simulate
posC      = find(strcmp(model.labely,'$c_t$'));
posW      = find(strcmp(model.labely,'$w_t$'));
pos_muz   = find(strcmp(model.labelx,'$\mu _{z,t}$'));
pos_n     = find(strcmp(model.labelx,'$n_{t}$'));
pos_d     = find(strcmp(model.labelx,'$d_{t}$'));

simC     = NaN(numSim,pathLength);
simZ     = NaN(numSim,pathLength);
simD     = NaN(numSim,pathLength);
simN     = NaN(numSim,pathLength);
simW     = NaN(numSim,pathLength);
simWstar = NaN(numSim,pathLength);
KAPAw    = model.params.KAPAw;
Wss      = exp(model.gSS(posW));

for i = 1:numSim
    %shocks = randn(size(model.eta,2),pathLength);
    outSim = simProjection(model,shocksMat(:,:,i));
    ySim   = outSim.y_cu;
    xSim   = outSim.xAll_cu(1:model.nx,:);
    if model.appOrder > 1
        xSim(1:model.nx1,:)   = xSim(1:model.nx1,:) + outSim.xAll_cu(model.nx+1:model.nx+model.nx1,:);
    end
    if model.appOrder > 2
        xSim(1:model.nx1,:)   = xSim(1:model.nx1,:) + outSim.xAll_cu(model.nx+model.nx1+1,model.nx+2*model.nx1);
    end
    
    muz       = model.hSS(pos_muz)   + xSim(pos_muz,:);
    n         = model.hSS(pos_n)     + xSim(pos_n,:);
    d         = model.hSS(pos_d)     + xSim(pos_d,:);
    simZ(i,:) = cumprod(exp(muz),2); 
    simN(i,:) = exp(n);
    simD(i,:) = exp(d); 
    simC(i,:) = exp(ySim(posC,:)).*simZ(i,:);
    simW(i,:) = exp(ySim(posW,:)).*simZ(i,:);
    simWstar(i,:) = (simW(i,:)-KAPAw*Wss*simZ(i,:))./(1-KAPAw);
end

bettap = params.BETTA.^(0:1:pathLength-1);

setupTP.N0  = params.U0*sum(repmat(bettap(1,1:end-1),numSim,1)...
              .*simZ(:,2:end).^((1-params.CHI)*(1-params.CHI0)),2);
Css         = exp(model.gSS(posC));
setupTP.Nc  = 1/(1-params.CHI)*sum(repmat(bettap(1,1:end-1),numSim,1).*simD(:,2:end)...
              .*((simC(:,2:end)-params.B*simC(:,1:end-1))./((Css*simZ(:,2:end)).^params.CHI0)).^(1-params.CHI),2);

setupTP.Nd  = params.U0d*sum(repmat(bettap(1,1:end-1),numSim,1).*simD(:,2:end)...
              .*simZ(:,2:end).^((1-params.CHI)*(1-params.CHI0)),2);

PHIzero     = model.params.PHIzero;
setupTP.Nl  = PHIzero^(params.PHI)/(1-1/params.PHI)*...
              sum(repmat(bettap(1,1:end-1),numSim,1).*simD(:,2:end).*simN(:,2:end).^params.PHI...
             .*simZ(:,2:end).^(params.PHI*(1-params.CHI)*(1-params.CHI0))...
             .*(simC(:,2:end)-params.B*simC(:,1:end-1)).^(params.CHI*(params.PHI-1))...
             .*(Css.*simZ(:,2:end)).^(params.CHI0*(params.PHI-1)*(1-params.CHI))...
             ./simWstar(:,2:end).^(params.PHI-1),2);
setupTP.c_cu           = Css;
setupTP.w_cu           = Wss;
setupTP.z_cu           = 1;
setupTP.n_cu           = 1;
setupTP.d_cu           = 1;
setupTP.params         = params;
setupTP.params.PHIzero = PHIzero;
setupTP.Css            = Css;
% Q = timingPremiumObjectFunc(0.5,setupTP);

options = optimset('Display','iter','MaxFunEvals',2000,'MaxIter',2000,'TolFun',1D-4,'TolX',1D-3);
paiOpt = fminsearch(@timingPremiumObjectFunc,pai0,options,setupTP);
%paiOpt    = fminunc(@timingPremiumObjectFunc,pai0,options,setupTP);

end